!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
program main
!==================================================================================================================================|

   use parameters
   use block_define
   use interfaces

   implicit none

   integer(4) :: ns, i
   real(8) :: t, t_start, t_stop, t_wall

   character(4) :: nt_char
   real(8) :: vx, vy, vz, b2, v2, pr, c2, s2, ca2, cfx

   type(block), pointer :: base

!----------------------------------------------------------------------------------------------------------------------------------|
   call cpu_time(t_start)

!----------------------------------------------------------------------------------------------------------------------------------|
   write (*, *) 
   write (*, '(a)') '|---------------|'
   write (*, '(a)') ' Program starts:'
   write (*, '(a)') '|---------------|'
   write (*, *)

!----------------------------------------------------------------------------------------------------------------------------------|
!   open (unit_in, file = 'input')
!   read (unit_in, '(a)') directory
!   read (unit_in, '(a)') parameters_filename
!   close (unit_in)

   call getcwd(directory)
   directory = trim(directory) // '/'

   call read_parameters
   call system("mkdir -p " // trim(datamap_path))

   open(unit_datamap_manage_module, file = 'data_manage_precision.include')
   write(unit_datamap_manage_module, '(a, i)') 'integer(4), parameter :: output_real8_flag = ', output_real8_flag
   close(unit_datamap_manage_module)
   open(unit_datamap_manage_module, file = 'data_manage_parameter.include')
   write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(parameters_filename)),                              &
      ') :: parameters_filename = ''', trim(parameters_filename), ''''
   write(unit_datamap_manage_module, '(a, i, a, a, a)') 'character(', len(trim(datamap_path)),                                     &
      ') :: datamap_path = ''', trim(datamap_path), ''''
   close(unit_datamap_manage_module)

!----------------------------------------------------------------------------------------------------------------------------------|
   t = 0
   t_wall = 0
   ns = 1

!----------------------------------------------------------------------------------------------------------------------------------|
   allocate (base)
   allocate (base%x(nx_total) ,base%ro(nx_total), base%mx(nx_total), base%my(nx_total), base%mz(nx_total),                         &
      base%bx(nx_total), base%by(nx_total), base%bz(nx_total), base%en(nx_total), base%child(ini_blocks),                          &
      base%flux_bnd(n_neighs, sys), base%neigh(n_neighs, n_levels))

   call model_initial (base)

!----------------------------------------------------------------------------------------------------------------------------------|
   open (unit_par, file = trim(datamap_path) // '/params.mapout')
   write (unit_par, '(2x,a)') trim(directory)
   write (unit_par, '(i5,t64,a)') ng, 'ng'
   write (unit_par, '(i5,t64,a)') nx_no_bnd, 'nx_no_bnd'
   write (unit_par, '(i5,t64,a)') nx_block, 'nx'
   write (unit_par, '(i5,t64,a)') nx_total, 'nx_total'
   write (unit_par, '(i5,t64,a)') ini_blocks, 'ini_blocks'
   write (unit_par, '(i5,t64,a)') ratio, 'ratio'
   write (unit_par, '(i5,t64,a)') n_levels, 'n_levels'
   write (unit_par, '(f10.5,t64,a)') gm, 'gm'
   write (unit_par, '(f10.5,t64,a)') rangex, 'rangex'
   write (unit_par, '(f10.5,t64,a)') zero_posx, 'zero_posx'
   write (unit_par, '(f10.5,t64,a)') cfl, 'cfl'
   close (unit_par)

!----------------------------------------------------------------------------------------------------------------------------------|
   do while (.true.)

      if (t .ge. nt * dt_output) then

         nt = nt + 1

         write(*,"(a5, i6, a6, i4, a11, i6, a5, e12.4, a6, e12.4)")                                                                &
            ' ns= ', ns, ', nt= ', nt, ', blocks = ', n_blocks, ', t= ', t, ', dt= ', min_dt

         write(nt_char,'(i4.4)') nt

         open (unit_nt, file = trim(datamap_path) // '/nt.mapout', form = 'unformatted')
         write (unit_nt) nt
         close (unit_nt)

         open (unit_t, file = trim(datamap_path) // '/info.t' // nt_char // '.mapout', form = 'unformatted')
         write (unit_t) t
         write (unit_t) min_dt
         close (unit_t)

         open(unit_datamap, file = trim(datamap_path) // '/data.t' // nt_char // '.mapout', form = 'unformatted')
         write (unit_datamap) n_blocks

         call block_output (base)
         close(unit_datamap)

         if (t .ge. t_end) exit

      endif

      call amr (base)

!      comment these lines for periodic boundary (the initial procedure should be change)
      call bnd_free_x (base%ro, ng, nx_total, 'b')
      call bnd_free_x (base%mx, ng, nx_total, 'b')
      call bnd_free_x (base%my, ng, nx_total, 'b')
      call bnd_free_x (base%mz, ng, nx_total, 'b')
      call bnd_free_x (base%bx, ng, nx_total, 'b')
      call bnd_free_x (base%by, ng, nx_total, 'b')
      call bnd_free_x (base%bz, ng, nx_total, 'b')
      call bnd_free_x (base%en, ng, nx_total, 'b')

      t = t + min_dt

      max_speed = 0.d0
      do i = ng + 1, nx_total - ng
         vx = base%mx(i) / base%ro(i)
         vy = base%my(i) / base%ro(i)
         vz = base%mz(i) / base%ro(i)
         b2 = base%bx(i) * base%bx(i) + base%by(i) * base%by(i) + base%bz(i) * base%bz(i)
         v2 = vx * vx + vy * vy + vz * vz
         pr = (base%en(i) - v2 * base%ro(i) / 2.0d0 - b2 / 2.0d0) * gmm1
         c2 = gm * pr / base%ro(i)
         s2 = c2 + b2 / base%ro(i)
         ca2 = base%bx(i) * base%bx(i) / base%ro(i)
         cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / 2.0d0)
         max_speed = max(max_speed, cfx + abs(vx))
      enddo

      base_dt = cfl * base%dx / max_speed

      ns = ns + 1

   enddo

   deallocate (base)

   call cpu_time (t_stop)

   t_wall = t_stop - t_start

   write (*,*)
   write (*,*) ' Total blocks: ', id_index
   write (*,"(a14,i3,a1,i3,a1,f6.2,a1)") " Elapsed time:", int(t_wall/3600.), 'h', int(mod(t_wall,3600.)/60), 'm', mod(t_wall,60.), 's'
   write (*,*)

   stop
end program main
